Write a custom CUDA kernel to optimize the fused `RMSNorm + SiLU` operator.

Formula:
1. RMSNorm: x_norm = x * w * rsqrt(mean(x^2) + eps)
2. SiLU: output = x_norm / (1 + exp(-x_norm))

Problem Analysis:
1. Memory Bound: A naive implementation executes RMSNorm (read input, calculate stats, write output) followed by SiLU (read output, compute, write final). This involves redundant memory round-trips.
2. Architecture Fit: Modern LLMs (Llama) use hidden dimensions like 4096. A single row (4096 * 4 bytes = 16KB) fits entirely in the Shared Memory of a standard CUDA block (usually 48KB+), enabling a "Cache-Once" strategy.

Optimization Strategy: Shared Memory Cached Fusion

1. Block-per-Row: Launch one CUDA block for each row (token) of the input batch.

2. Shared Memory Caching:
   - Cooperative Load: Threads cooperatively load the entire row from Global Memory into Shared Memory using `float4` vectorized loads.
   - This ensures the input is read from HBM only once.

3. Two-Pass Algorithm (in Shared Memory):
   - Pass 1 (Reduction): Compute `sum(x^2)` using parallel reduction in Shared Memory. Calculate `inv_rms = rsqrt(sum/N + eps)`.
   - Pass 2 (Apply): Read value from Shared Memory, apply normalization, multiply by weight gamma, apply SiLU activation, and store back to Global Memory.

4. Vectorization: Use `float4` for all Global Memory accesses and `float` logic for computation. Ideally process multiple elements per thread to hide latency.

Here's an example to show you the syntax of inline embedding custom CUDA operators in torch: The example given architecture is:   
  
```python
import torch
import torch.nn as nn

BATCH_SIZE = 4096
HIDDEN_DIM = 4096
SHAPE = (BATCH_SIZE, HIDDEN_DIM)

EPS_VALUE = 1e-5

class RMSNorm(nn.Module):
    """
    RMSNorm 以兼容旧版 PyTorch。
    """
    def __init__(self, hidden_size, eps=1e-5):
        super(RMSNorm, self).__init__()
        self.weight = nn.Parameter(torch.ones(hidden_size))
        self.variance_epsilon = eps

    def forward(self, x):
        input_dtype = x.dtype
        x = x.to(torch.float32)
        
        variance = x.pow(2).mean(-1, keepdim=True)
        
        hidden_states = x * torch.rsqrt(variance + self.variance_epsilon)
        
        return (self.weight * hidden_states).to(input_dtype)

class RMSNormSiLU(nn.Module):
    """
    Fused RMSNorm + SiLU
    """
    def __init__(self, dim, eps=1e-5):
        super(RMSNormSiLU, self).__init__()
        self.rms = RMSNorm(dim, eps=eps)
        self.silu = nn.SiLU()

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        # RMSNormKernel -> SiluKernel
        x = self.rms(x)
        return self.silu(x)

class Model(nn.Module):
    def __init__(self, dim, eps=1e-5):
        super(Model, self).__init__()
        self.fused_op = RMSNormSiLU(dim, eps)
    
    def forward(self, x):
        return self.fused_op(x)

def get_inputs():
    input_tensor = torch.randn(SHAPE, dtype=torch.float32)
    return [input_tensor.contiguous()]

def get_init_inputs():
    return [HIDDEN_DIM, EPS_VALUE]